function [beta,se,p,r2,F,control_mean,N,C,T] = fereg(y,treat,X,stations,tail);
    %%Estimation Sample%%
    [N,k] = size(X);
    T = sum(treat);
    C = sum(1-treat);
    
    control_mean = sum((1-treat).*y)/sum((1-treat));
    
    Xcov = [treat X];
    [beta,bint,e,rint,stats] = regress(y,Xcov);
      
    dep_matrix = depmat(N,stations);
        
    omega = (e*e').*dep_matrix;
    V = inv(Xcov'*Xcov)*(Xcov'*omega*Xcov)*inv(Xcov'*Xcov);
    se = sqrt(diag(V));
    
    sst = sum((y-mean(y)).^2);
    ssm = sum((Xcov*beta-mean(y)).^2);
    sse = sum(e.^2);
   
    r2 = stats(1);
    F = (ssm/k)/(sse/(57-k-1));
    
    p = zeros(k+1,1);
    
    for i = 1:k+1
        if tail == 1
            p(i) = 1 - cdf('norm',abs(beta(i)/se(i)),0,1);
        elseif tail == 2 
            p(i) = 2*(1 - cdf('norm',abs(beta(i)/se(i)),0,1));
        end;
    end;

    
    
        
    
    
    
    
    
    